#pragma once

// 数值计算的精度最好取 double，不然精度不够
using Real = double;

#include <cmath>
#include <climits>
#include <cfloat>
#include <cassert>
#include <initializer_list>

namespace LinearAlgebra
{
	// 预先声明
	class Vector;
	class Matrix;

	/**
	 * @file Matrix.h
	 * @author xingyifan
	 * @date 2022-11-28 21:51
	 *
	 * @description: 向量类型
	 */

	class Vector
	{
	protected:
		const int m_dim;
		Real *m_value;

	public:
		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 21:54
		 *
		 * @description: 构造函数
		 */

		Vector(int D, Real t = 0);			   // 初始化构造函数
		Vector(std::initializer_list<Real> l); // 可以通过列表初始化
		Vector(const Vector &rhs);			   // 拷贝构造

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-30 09:25
		 *
		 * @description: 向量信息
		 */

		int size() const; // 获得向量长度

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 21:57
		 *
		 * @description: 操作符重载
		 */

		Real &operator[](int i);	  // 返回引用
		Real operator[](int i) const; // 返回值

		// 向量加减
		Vector operator+(const Vector &rhs) const;
		Vector operator-(const Vector &rhs) const;

		// 常数乘除法
		Vector operator*(Real t) const;
		Vector operator/(Real t) const;

		// 赋值
		void operator=(const Vector &rhs);

		// 取反操作
		Vector operator-() const;

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 22:12
		 *
		 * @description: 基本向量运算
		 */

		void sort();				 // 冒泡排序
		Vector reverse();			 // 反转向量
		Real sum() const;			 // 元素求和
		Vector abs() const;			 // 取绝对值
		Real norm(int nt = 2) const; // 向量范数
		Real Min() const;			 // 取最小元
		Real Max() const;			 // 取最大元
		int maxIndex() const;		 // 取具有最大模的分量
		Matrix toMatrix() const;	 // 转化成矩阵
		Matrix transpose() const;	 // 转置成矩阵

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-30 09:27
		 *
		 * @description: 释放资源
		 */

		~Vector();
	};

	/**
	 * @file Matrix.h
	 * @author xingyifan
	 * @date 2022-11-28 22:29
	 *
	 * @description: 矩阵类型
	 */

	class Matrix
	{
	protected:
		const int m_row, m_col;
		Real **m_value;

	public:
		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 22:30
		 *
		 * @description: 构造函数
		 */

		Matrix(int Row, int Col, Real t = 0);					 // 初始化构造函数
		Matrix(int Row, int Col, std::initializer_list<Real> l); // 通过列表初始化
		Matrix(const Matrix &rhs);								 // 拷贝构造

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-30 09:52
		 *
		 * @description: 矩阵信息
		 */

		int row() const; // 获得行数
		int col() const; // 获得列数

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 22:37
		 *
		 * @description: 操作符重载
		 */

		Vector operator()(int j) const;		 // 返回第 j 列
		Real operator()(int i, int j) const; // 返回 i,j 元素
		Real &operator()(int i, int j);		 // 返回 i,j 元素引用

		// 矩阵加减法，矩阵乘向量
		Matrix operator+(const Matrix &rhs) const;
		Matrix operator-(const Matrix &rhs) const;
		Matrix operator*(Real t) const;
		Matrix operator/(Real t) const;
		Vector operator*(const Vector &v) const;
		Matrix operator*(const Matrix &rhs) const;

		void operator=(const Matrix &rhs); // 直接赋值
		Matrix operator-() const;		   // 取反操作

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-28 23:05
		 *
		 * @description: 元素操作
		 */

		void insert(int j, Vector &v); // 在 j 列插入列向量 v
		Real norm(int nt = 0) const;   // 矩阵范数
		Real norm(char c) const;	   // e 范数 / f 范数
		Matrix upper() const;		   // 获得上三角
		Matrix lower() const;		   // 获得下三角
		Matrix transpose() const;	   // 转置

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-30 10:34
		 *
		 * @description: 矩阵算法，不改变矩阵
		 */

		Matrix inverse() const;				// 矩阵求逆
		Matrix LU() const;					// 标准 LU 分解
		Matrix LU(int perm[]) const;		// 列主元 LU 分解，传入下标数组，返回置换后的数组
		Matrix Cholesky() const;			// Cholesky LL'
		Matrix CholeskyD() const;			// 优化的 Cholesky LDL'
		Matrix QR(Vector &d) const;			// QR 分解
		Matrix Hessenberg(Matrix &Q) const; // 上 Hessenberg 化，在 Q 中累计正交变换

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-12-01 15:32
		 *
		 * @description: 矩阵迭代法，将对矩阵本身进行迭代
		 */

		void HAH(int startOfRow, int endOfRow, Vector &v); // 计算指定范围的 Householder 变换，作用结果为 HAH
		void doubleQR(int l = 0, int m = 0);			   // 双重步 QR 迭代，l,m 是 H22 前面和后面分块的阶数

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-30 09:40
		 *
		 * @description: 释放资源
		 */

		~Matrix();
	};

	/**
	 * @file Matrix.h
	 * @author xingyifan
	 * @date 2022-11-29 09:39
	 *
	 * @description: 动态数组，用于稀疏矩阵
	 */

	template <class T>
	class DynamicArray
	{
	public:
		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 09:39
		 *
		 * @description: 内置节点类型
		 */

		struct Node
		{
			T value;
			Node *prev;
			Node *next;

			Node(const T &value = T(), Node *prev = nullptr, Node *next = nullptr) : prev(prev), next(next), value(value) {}
		};

	protected:
		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 09:40
		 *
		 * @description: 保有节点头部和数组长度
		 */

		int m_size;
		Node *m_head;
		Node *m_tail;

	public:
		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 09:41
		 *
		 * @description: 构造函数
		 */

		DynamicArray();							  // 直接构造
		DynamicArray(std::initializer_list<T> l); // 初始化构造
		DynamicArray(const DynamicArray &arr);	  // 拷贝构造

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 10:26
		 *
		 * @description: 获得数组信息
		 */

		int size() const { return m_size; }
		bool empty() const { return (m_size == 0) ? true : false; }

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 10:04
		 *
		 * @description: 操作符重载
		 */

		T &operator[](int i);
		T operator[](int i) const;

		// 拷贝赋值
		void operator=(const DynamicArray &arr);

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 10:09
		 *
		 * @description: 元素操作
		 */

		T front() const;				   // 获得头部元素
		T back() const;					   // 获得尾部元素
		Node *node(int i) const;		   // 获得指定的节点
		void insert(int i, T value);	   // 在 i 位置插入元素
		Node *insert(Node *node, T value); // 在指定节点处插入元素并返回新节点
		void push_back(T value);		   // 最后加入元素
		void push_front(T value);		   // 最前面加入元素
		T pop_back();					   // 弹出最后一个元素
		T pop_front();					   // 弹出第一个元素
		int find(T value) const;		   // 找到元素
		bool kill(T value);				   // 删除指定元素
		T remove(int i);				   // 移除指定索引处的元素并返回
		void remove(Node *node);		   // 移除指定节点的元素
		void clear();					   // 删除数组元素

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 10:27
		 *
		 * @description: 释放内存
		 */

		~DynamicArray();
	};

	/**
	 * @file Matrix.h
	 * @author xingyifan
	 * @date 2022-11-29 14:09
	 *
	 * @description: 稀疏矩阵类型
	 */

	class SpMatrix
	{
	protected:
		const int m_row, m_col;		  // 尺寸
		DynamicArray<int> m_IA, m_JA; // 存放索引
		DynamicArray<Real> m_AA;	  // 对角元放在每行首位

	public:
		friend class LinearSolver; // 定义友元，方便访问

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 14:10
		 *
		 * @description: 构造函数
		 */

		SpMatrix(int Row, int Col);	   // 直接构造
		SpMatrix(const SpMatrix &spm); // 拷贝构造
		SpMatrix(const Matrix &m);	   // 用普通矩阵构造

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 15:45
		 *
		 * @description: 结构信息
		 */

		int row() const;				   // 获得行数
		int col() const;				   // 获得列数
		int explicitNum() const;		   // 获得显式表达元素的个数
		Real diag(int i) const;			   // 获取对角元
		int startOfRow(int i) const;	   // 获取指定行的开头
		int indexOfCol(int j) const;	   // 通过全局指标获取元素所在的列位置
		Real global(int j) const;		   // 通过全局指标获取元素值
		int globalCol(int i, int j) const; // 获取(i,j)元素的全局列标

		/**
		 * @file Matrix.h
		 * @author xingyifan
		 * @date 2022-11-29 15:51
		 *
		 * @description: 插入和修改元素
		 */

		Real operator()(int i, int j) const; // 获得存在的元素，不存在的元素认为是零
		Real &operator()(int i, int j);		 // 获得元素引用，不存在的元素就创建元素

		SpMatrix compress() const;	// 压缩矩阵的零元素，返回压缩矩阵
		SpMatrix transpose() const; // 转置矩阵
		Matrix toMatrix() const;	// 转化为一般矩阵

		/**
		 * @file Matrix.hpp
		 * @author xingyifan
		 * @date 2022-11-29 17:14
		 *
		 * @description: 操作符重载
		 */

		// 加减法
		SpMatrix operator+(const SpMatrix &spm) const;
		SpMatrix operator-(const SpMatrix &spm) const;

		// 乘法运算
		SpMatrix operator*(Real k) const;
		SpMatrix operator/(Real k) const;

		// 稀疏矩阵乘稀疏矩阵
		SpMatrix operator*(const SpMatrix &spm) const;

		// 拷贝赋值
		void operator=(const SpMatrix &spm);
	};

	/**
	 * @file Matrix.h
	 * @author xingyifan
	 * @date 2022-11-30 21:36
	 *
	 * @description: 基本向量操作和矩阵操作
	 */

	Real dot(const Vector &vec1, const Vector &vec2); // 向量点积
	Real Householder(Vector &x, Vector &v);			  // Householder 变换
	Vector Givens(Vector &v, int i, int j);			  // Givens 变换
	Matrix Hilbert(int n);							  // 希尔伯特矩阵
};